Stabilisation of the lattice-Boltzmann method 
using the Ehrenfests' coarse-graining 
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The lattice-Boltzmann method (LBM) and its variants have emerged as promising, computation- 
ally efficient and increasingly popular numerical methods for modelling complex fluid flow. However, 
it is acknowledged that the method can demonstrate numerical instabilities, e.g., in the vicinity of 
shocks. We propose a simple and novel technique to stabilise the lattice-Boltzmann method by 
monitoring the difference between microscopic and macroscopic entropy. Populations are returned 
to their equilibrium states if a threshold value is exceeded. We coin the name Ehrenfests' steps 
for this procedure in homage to the vehicle that we use to introduce the procedure, namely, the 
Ehrenfests' idea of coarse-graining. 
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I. INTRODUCTION 

The lattice-Boltzmann method (LBM) provides an al- 
ternative to the orthodox approach to computational 
fluid dynamics, in which the starting point is always 
a discretization of the Navier-Stokes equations. The 
method, which is fundamentally based on Boltzmann's 
kinetic transport equation, instead describes a fluid by 
a number of interacting populations of particles moving 
and colliding on a fixed lattice. With the advent of the 
introduction of a diagonal collision integral with single- 
time relaxation to local equilibrium, the method becomes 
tantalisingly simple and efficient. 

In recent years, the LBM has enjoyed much applied 
success modelling various flows of genuine engineering in- 
terest (see e.g., l] and the references therein). However, 
when populations are far from equilibrium, such as is the 
case in the vicinity of shocks, the LBM exhibits numer- 
ical instability. Often, numerical instability in the LBM 
is attributed to the absence of a positivity constraint 
on the populations. Recent development of the entropic 
LBM (ELBM) 0, H are attem P ts to improve stability 
properties through compliance with a discrete entropy H- 
theorem. Although such techniques should, conceivably, 
reduce the magnitude of spurious oscillations, one should 
perhaps not be disappointed to find that such numerical 
instabilities are not removed entirely. 

In this paper we suggest an alternative and versatile 
approach. The idea is simply stated: we propose a LBM 
in which the difference between microscopic and macro- 
scopic entropy is monitored in the simulation, and popu- 
lations are returned to their equilibrium states if a thresh- 
old value is exceeded. This technique is appealing be- 
cause dissipation is introduced in a controlled, targeted 
and illiberal manner. Furthermore, we stress that equi- 
libration itself will leave macroscopic entropy completely 



unchanged. 

We coin the name Ehrenfests ' steps for this procedure 
because we feel that the technique is most clearly under- 
stood if introduced using the Ehrenfests' coarse-graining 
idea 

This report is organised as follows. In Section ITT1 the 
LBM is recalled. In Section ILTTI we recall the basic theory 
required to discuss the Ehrenfests' coarse-graining and 
introduce Ehrenfests' steps. The reader is directed to Q 
for a more comprehensive review of coarse-graining. Fi- 
nally, we present the results of a shock tube numerical 
experiment which compares various LBMs. 



II. LATTICE-BOLTZMANN METHOD 

The Boltzmann kinetic transport equation is the fol- 
lowing time evolution equation for one-particle distribu- 
tion functions / = f(x,v,t): 



d t f 



V/ = Q. 



The collision integral, Q, describes the interactions of the 
populations /. The lattice-Boltzmann approach drasti- 
cally simplifies this model by stipulating that popula- 
tions can only move with a finite number of velocities 
{vi, ■ ■ -,v n }: 



dtfi 



V/i = Qi, i = 1, . 



(1) 
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where /, is the one-particle distribution function asso- 
ciated with motion in the ith direction. For example, 
in one-dimension, one might consider three velocities 
{— c, 0, c}, for some c ^ 0. 

The model is further simplified by specifying the colli- 
sion integral as the Bhatnager-Gross-Krook (BGK) op- 
erator Q: Qi = (/j 6q — fi)/r. We make this assumption 
for the remainder of the paper. Now, (JTJ describes free- 
flight dynamics together with relaxation to local equilib- 
ria, / oq , in time proportional to r. 
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The discrete velocities and local equilibrium states 
can all be mindfully chosen so that the Navier-Stokes 
equations are recovered by the lattice-Boltzmann equa- 
tion subject to ceratin conservation laws, in the large 
time-scale, t, limit via a Chapman-Enskog procedure 
Here, the macroscopic fluid density p and momentum pu 
are the zeroth- and first-order hydrodynamic moments 
of the populations respectively. The rate of dissipation 
introduced by the model is proportional to r. 

The local equilibrium states can be found by maximis- 
ing a local entropy functional subject to the constraints 
of conservation of mass and momentum. 

There are numerous ways to discretizc and obtain a 
numerical method. We prefer the following, second-order 
accurate in time, lattice-based LBGK scheme 0: 

f i (x + v i At,t + At) = (l-0)Mx,t)+0Mx,t), (2) 

with fi — 2/j Cq — /j. Here, the discrete velocities are as- 
sociated with an underlying spatial lattice L and x G C. 
Populations live on this lattice, propagate to neighbour- 
ing lattice sites with their corresponding discrete veloc- 
ities and are updated via J5J. The parameter (3 G (0, 1] 
controls the viscosity in the model, with [3 = 1 corre- 
sponding to the zero-viscosity limit. 

III. THE EHRENFESTS' COARSE-GRAINING 

To introduce the Ehrenfests' coarse-graining idea Q 
we use a formal kinetic equation 

%=J(f), 0) 

together with a strictly concave entropy functional S. 
We tacitly assume that entropy does not decrease with 
time: dS(f)/dt > 0. For our purposes, we are always 
interested in the example of the lattice-Boltzmann equa- 
tion Q with associated entropy functional which defines 
the local equilibrium states. 

The Ehrenfests' idea was to supplement the mechan- 
ical motion from @ with periodic averaging in cells to 
produce piecewise constant, or coarse-grained, densities. 
This operation necessitates entropy production. 

We wish to allow a generalisation of the Ehrenfests' 
coarse-graining whereby averaging in cells is replaced 
with some other partial equilibration procedure. Specif- 
ically, we assume we have some linear operator m which 
transforms a microscopic description of the system / into 
a macroscopic description M = m(/). For our example, 
the macroscopic description is that provided by the usual 
hydrodynamic moments. Now, given a macroscopic de- 
scription, M, we consider the solution flj of the optimi- 
sation problem 

argmax{s(/) : m(/) = ikf}. (4) 

Averaging in cells is a particular example of this en- 
tropy maximisation problem for the Boltzmann-Gibbs- 
Shannon (BGS) entropy functional S(f) = — J /log(/), 



where the integration is taken over the whole of phase 
space. 

We will refer to f^j as the quasi-equilibrium distri- 
bution. The quasi-equilibrium manifold Q is the set 
of quasi-equilibrium distributions parameterised by the 
macroscopic variables M. 



A. The Ehrenfests' chain and entropic involution 

Entropy maximisation leads naturally to an evolution 
equation for the macroscopic description. The so-called 
quasi-equilibrium approximation to © is an equation for 
the evolution of M: 

^=rn(J(f* M )). (5) 

For our example, the quasi-equilibrium distributions are 
precisely the local equilibria / 4 eq and 10 coincides with 
the compressible Euler equations Q . 

Now, one can envisage constructing various coarse- 
graining chains which provide stepwise approximation to 
the macroscopic equations (0. 

Let 0t be the phase flow for the kinetic equation 0. 
Let r be a fixed coarse- graining time and suppose we have 
an initial quasi-equilibrium distribution / . The Ehren- 
fests' chain is the sequence of quasi-equilibrium distribu- 
tions / ,/i, . . ., where ft := f^e^ft-t))' 

Entropy increases in the Ehrenfests' chain. The en- 
tropy gain in a link in the chain is made up from two 
parts: the entropy gain from the mechanical motion 
(from fi to Q T (fi)) and the gain from the equilibration 
(from T (/i) to fi+i). Consequently, conservative sys- 
tems become dissipative and dissipative systems more 
so. The gain in macroscopic entropy in a particular link 
is given by the expression S{f* n[h+i) ) - S{f* nUi) ). Note 
that there is zero gain in macroscopic entropy from the 
equilibration part. 

The Ehrenfests' chain provides a stepwise approxima- 
tion to a solution of some coarse-grained macroscopic 
equations via Mi := m{fi), i = 1,2,.... For our exam- 
ple, these equations are the compressible Navier-Stokes 
equations However, this chain is computationally 
prohibitive because the rate of introduced dissipation is 
proportional to r. The Ehrenfests' chain corresponds to 
the following LBM: 

fi(x + Vi At,t + At) = {l-p)fi(x,t)+f3fP{x,t), 

which we recognise as a forward Euler discretization 
of p. 

Another possibility is to construct a chain as follows. 
As already noted, the dissipative term introduced by the 
Ehrenfests' chain depends linearly on r. Therefore, there 
is symmetry between forward and backward motion in 
time starting from any quasi-equilibrium initial condi- 
tion. It is precisely this principle that enables one to 
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construct chains with zero macroscopic entropy produc- 
tion. Each subsequent link in the chain is constructed 
using entropic involution. 

To make such a chain useful, the user is at liberty to 
add a required amount of dissipation by shifting the invo- 
luted point in the direction of the quasi-equilibrium state, 
with some entropy increase. Of course, this shift leaves 
macroscopic entropy unchanged. This chain corresponds 
to the ELBM that we have eluded to in the introduction: 



f l (x + v l At,t + At) = {l-0)fi(x,t)+0fi, a (x,t), (6) 

with f i>a = (1 - a)fi + af° q . The number a = a(fi) 
is chosen so that a constant entropy condition is satis- 
fied. This, in itself, provides a positivity constraint on 
the populations. 

If the entropic involution operator is partially lin- 
earised so that subsequent links are constructed with 
the use of reflections in the quasi-equilibrium manifold, 
rather than inversions, then the corresponding LBM is 
precisely LBGK ©. 



B. Ehrenfests' steps 

We can now state the main idea of the paper. We 
propose to create a new chain. This chain begins and 
proceeds, for the bulk of time, as the either the entropic 
involution chain or its linearised version, as described 
in the previous subsection. However, the difference be- 
tween microscopic and macroscopic entropy is monitored 
throughout the simulation, by which we mean the quan- 
tities 



AS, := £(/, 



m(/»; 



S(fi). 



A threshold value is set and an alarm is triggered if ex- 
ceeded. The alarm simply signals that a link from the 
Ehrenfests' chain — an Ehrenfests' step — be used in 
place of a regular link of the primary chain at this point. 
Links in the chain which are intolerably far from their 
quasi-equilibrium states are merely returned to their 
equilibrium. The result is a chain which provides ad- 
ditional dissipation only where it is anticipated to be re- 
quired. 

Coupling the entropic involution chain with Ehren- 
fests' steps will, in general, no longer constitute a chain 
with zero macroscopic entropy production. Indeed, in a 
chain of length N + 1 the total macroscopic entropy gain 
is given by the expression £ ieJ [S '(fm(f i+1 )) ~ S r (/m(/,))] ' 
where I C {1, . . . , N} is the set of indices corresponding 
to the Ehrenfests' steps. 



IV. NUMERICAL EXPERIMENT 

To conclude this report, we perform a one-dimensional 
numerical experiment to demonstrate the performance 



of the proposed LBM corresponding to a chain with 
Ehrenfests' steps. The implementation is as follows. If 
A Si > S, for some given tolerance 5 > then we accept 
an Ehrenfests' step. The resulting LBM is as follows: 



fi(x + Vi£t,t + At) 



(l-/3)/i + A/T. 
(1 -/?)/; +/?/;, 



AS, > 5, 

otherwise. 
(7) 

We have selected LBGK as the primary chain and we 
will henceforth refer to this LBM as LBGK-ES. 

The method LBGK-ES, as described, is a second-order 
accurate in time scheme with first-order degradation in 
regions where Ehrenfests' steps are employed. However, 
as is often found, when the thickness of such regions is 
of the order of the lattice spacing, the method remains 
second-order everywhere. 

For contrast we are interested in comparing LBGK-ES 
with LBGK (JTJ and ELBM © (as described in 0). In 
all of our simulations we select the three-velocity model 
mentioned in the introduction and a uniformly spaced 
lattice. Further, we always choose c = 1 and At = 1. 
The coefficient (3, which controls the viscosity in the 
model, is fixed at /? = 1 — 10~ 9 , which is close to the 
zero viscosity limit. In each case the entropy is S = —H, 
with H = ft log(A/4) + h log(/ 2 ) + h log(/ 3 ) 0, where 
fit fa and fa denote the static, left-moving and right- 
moving populations respectively. For this entropy, the 
local equilibria are available analytically for the three- 
velocity model. They are given by the expressions 



/ 1 cq = 2p(2- v / T+3^)/3, 



J 2 cq = p((3u - 1) + 2 v / l + 3zt 2 )/6, 
(-(3u + 1) + 2Vl + 3u 2 )/6. 



J3 



For LBGK-ES we fix the tolerance, S, to either 5 = 10 3 
or 5 = 1(T 5 . 

As already mentioned, for ELBM, there is a parameter, 
a, which is chosen to satisfy a constant entropy condition. 
This involves finding the nontrivial root of the equation 



H(f + aQ) = H(f). 



(8) 



Inaccuracy in the solution of this equation can introduce 
artificial viscosity. To solve © numerically we employ 
a robust routine based on bisection. The root is solved 
to an accuracy of 10 -15 . Furthermore, we always ensure 
that the returned value of a does not lead to an entropy 
decrease. So that our results can be faithfully repro- 
duced, we stipulate that we will consider two possibilities 
if no nontrivial root of |JSJ exists. We either elect to se- 
lect a = 1 or we select the positivity bound a — a+ , with 
a + := muiQ^o \fi/Qi\- For fairness, since both of these 
procedures introduce diffusivity into ELBM, we consider 
analogous procedures for LBGK too, i.e., we elect to se- 
lect either a = 1 or a = a + if a populations is predicted 
to become negative. 
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FIG. 1: Density profile of the isothermal 1 : 2 shock tube 
simulation after 300 time steps using (a) LBGK @; (b) 
ELBM ©; (c) LBGK-ES Q with 8 = 10" 3 ; (d) LBGK- 
ES with S = 10" 5 . In this example, LBGK does not 
produce a negative population so the aforementioned regu- 
larisation procedures are redundant. Similarly, for ELBM in 
this example, the entropy estimate equation always has a non- 
trivial root. Sites where Ehrenfests' steps are employed are 
indicated by crosses. 



A. Shock tube results 

The one-dimensional shock tube for a compressible 
isothermal fluid is a standard benchmark test for hydro- 
dynamic codes. Our computational domain will be the 
interval [0, 800] and we discretize this interval with 801 
uniformly spaced lattice sites. We choose the initial den- 
sity ratio as 1 : 2. Initially, for x < 400 we set p = 1.0 
else we set p = 0.5. 

We observe that, of all the LBMs considered in the 
experiment, only the method which includes Ehrenfests' 
steps is capable of suppressing spurious post-shock oscil- 
lations (Fig. pi . 

The code used to produce the simulations in this sec- 
tion is freely available by making contact with the corre- 
sponding author. 

Conclusion 

Ehrenfests' steps introduce additional dissipation lo- 
cally, on the base of point-wise analysis of non- 
equilibrium entropy. They admit a huge variety of gen- 
eralisations: incomplete Ehrenfests' steps, partial invo- 
lution, etc. [6|. In order to preserve the second-order of 
LBM accuracy it is worthwhile to perform Ehrenfests' 
steps on only a small share of sites (the number of sites 
should be 0(Nh/£), where N is the number of sites, I is 
the macroscopic characteristic length and h is the lattice 
step) with highest ASi > 5. Numerical experiments show 
that even small shares of such steps drastically improve 
stability. More tests are presented in |To) . 
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FIG. 2: Density profile of the isothermal 1 : 2 shock tube simulation after 500 time steps using a) LBGK; b) LBGK-ES with 
S = 10 -5 ; c) LBGK with regularisation provided by selecting the positivity bound a = a+ if a population is predicted to 
become negative; d) LBGK with regularisation provided by selecting a = 1 if a population is predicted to become negative; e) 
ELBM with a — a+ if no nontrivial root of JSJ exists; f) ELBM with a = 1 if no nontrivial root of (|HJ exists. In this example, 
LBGK does not produce a negative population so the aforementioned regularisation procedures are redundant. Similarly, for 
ELBM in this example, the entropy estimate equation always has a nontrivial root. 



6 




FIG. 3: Density profile of the isothermal 1 : 10 shock tube simulation after 350 time steps using a) LBGK; b) LBGK-ES with 
S = 10 -5 ; c) LBGK with regularisation provided by selecting the positivity bound a = a+ if a population is predicted to 
become negative; d) LBGK with regularisation provided by selecting a = 1 if a population is predicted to become negative; e) 
ELBM with a = a+ if no nontrivial root of (SJ exists; f) ELBM with a = 1 if no nontrivial root of JH) exists. In this example, 
we acknowledge that the initial 1 : 10 density discontinuity produces a shock which is outside of the hydrodynamic regime of 
the three-velocity model that we are using. However, this experiment is useful for testing the computational stability of the 
various LBMs because populations are capable of becoming negative. 
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FIG. 4: Density profile of the isothermal 1 : 2 shock tube simulation after 300 time steps using LBGK-ES with a) 8 — 10 -2 ; 
b) 8 = 10~ 3 ; c) 8 = 10 -4 ; d) 8 = 10 -5 ; e) S = 10 -6 ; f ) 5 = 10~ 7 . We are using crosses to indicate where an Ehrenfests' step is 
being used in the simulation. As we expect, the profile becomes progressively more smoothed as 8 decreases. 
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FIG. 5: Density profile of the isothermal 1 : 2 shock tube simulation after 400 time steps using LBGK-ES with S = 10 -5 . 
Here, we are using the (k, <5)-rule. This rule specifies that only the k sites with the highest values of ASi > S are accepted. 
The simulation shows the result for a) k=0; b) k=l; c) k—2; d) k—4; e) k=8; f) k=16. We are using crosses to indicate 
where an Ehrenfests' step is being used in the simulation. We observe that just a few Ehrenfests' steps greatly reduce spurious 
oscillations. 
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FIG. 6: Profile of a 7^ 1 in the isothermal 1 : 2 shock tube simulation after 300 time steps using ELBM-ES with a) 
b) 8 = 10~ 3 ; c) 8 = 10" 4 ; d) 8 = 10" 5 ; e) 8 = 10" 6 ; f) 8 = 10" 7 . The gaps in the profiles occur precisely where 
have been removed so that the diminishing variation in a as 8 decreases can be more readily observed. 
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FIG. 7: Advection of a square profile after 1000 time steps using a) LBGK; b) ELBM and c) LBGK-ES with 8 = 10~ 3 . All 
with P = 0.999. 
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FIG. 8: Advection of a square profile after 1000 time steps using LBGK-ES (J3 = 0.999) with a) 8 = 10~ 2 ; b) 8 = 10~ 3 ; c) 
8 = 10~ 4 ; d) S = 10~ 5 ; e) 5 = 10~ 6 ; f) 5 = 10~ 7 . We are using crosses to indicate where an Ehrenfests' step is being used in 
the simulation. 
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FIG. 9: Advection of a square profile after 1000 time steps using LBGK-ES (f3 = 0.999) with 5 — 10~ 5 . Here, we are using 
the (k, 5)-rule. This rule specifies that only the k sites with the highest values of ASi > 5 axe accepted. The simulation shows 
the result for a) k=0; b) k=l; c) k=2; d) k=4; e) k=8; f) k=16. We are using crosses to indicate where an Ehrenfests' step is 
being used in the simulation. 



